%-----------------------------------------
% m_equivsolve.m
% Main computation for spatial equilibrium
% compute dynamic spatial equilibrium
% ----------------------------------------

% ************************
% *** Basic parameters ***
% ************************
psi=0.1; 
mu=0.15; 
gammaL=0.7; 
gammaH=0.1;
epsilon=8; 
eta=4; 

% *******************
% *** initial set ***
% *******************
Pop=1; 

R1_bar=0; 
R2_bar=Pop-R1_bar;
R_bar=[R1_bar;R2_bar]; 

L1_bar=0;
L2_bar=Pop-L1_bar;
L_bar=[L1_bar;L2_bar]; 

Lcom_bar=[0 0; 0 Pop]; 
H_bar=[0; 0.1]; 
% *******************
% *** solve model ***
% *******************
maxit=1000000; 
tol=1e-5; 
it=0; 
update=0.05; 

while it<maxit   
    w=(AA.^(1/gammaL)).*(Q_i.^(-gammaH/gammaL)).*(L_i.^(alpha/gammaL)); 
    lambda=Kmat.*(repmat(BB',2,1).^epsilon).*repmat(((R_i.^(beta*epsilon)).*(Q_i.^(-mu*epsilon)))',2,1).*repmat(w.^epsilon,1,2); 
    lambda=lambda./sum(sum(lambda)); 
    Lcom=(1-theta).*Lcom_bar+theta.*lambda; 
    I=sum(Lcom.*repmat(w,1,2),1)'; 
    R=(1-theta).*R_bar+theta.*sum(lambda,1)'; 
    L=(1-theta).*L_bar+theta.*sum(lambda,2); 
    R_err=abs(R-R_i); 
    L_err=abs(L-L_i); 

    % floor space prices 
    H=(1-psi).*H_bar+Q_i.^eta; 
    Q=(mu.*I+(gammaH/gammaL).*w.*L_i)./H;
    Q_err=abs(Q-Q_i); 
    
    if R_err<tol & L_err<tol & Q_err<tol 
        break
    else
    if rem(it,50)==1
        disp(['Largest Error is:',' ',num2str(it), ' ',num2str(max(R_err)), ' ',num2str(max(L_err))]);
    end
    R_i=(1-update)*R_i+update*R; 
    L_i=(1-update)*L_i+update*L;
    Q_i=(1-update)*Q_i+update*Q;
    it=it+1;
    end
end 


